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Random Matrix Theory (RMT) is a powerful statistical tool to model spectral fluctuations. This 
approach has also found fruitful application in Quantum Chromodynamics (QCD). Importantly, 
RMT provides very efficient means to separate different scales in the spectral fluctuations. We try 
to identify the equivalent of a Thouless energy in complete spectra of the QCD Dirac operator for 
staggered fermions from SU(2) lattice gauge theory for different lattice size and gauge couplings. 
We focus on the bulk of the spectrum. In disordered systems, the Thouless energy sets the universal 
' scale for which RMT applies. This relates to recent theoretical studies which suggest a strong 

' analogy between QCD and disordered systems. The wealth of data allows us to analyze several 

, statistical measures in the bulk of the spectrum with high quality. We find deviations which allows 

us to give an estimate for this universal scale. Other deviations than these are seen whose possible 
origin is discussed. Moreover, we work out higher order correlators as well, in particular three-point 
correlation functions. 
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> ■ I. INTRODUCTION 

cn ■ 

, It is now well established that Random Matrix Theory (RMT) accurately models spectral fluctuations in an abun- 

\^ ' dant variety of different systems, such as chaotic, disordered and many-body systems, see the review in Ref. 0|. In 

' recent years, RMT has in addition been successfully introduced into the study of certain aspects of Quantum Chro- 

00 ■ modynamics (QCD). The interest focuses on the spectral properties of the Euclidean Dirac operator. The eigenvalue 

[ equation under consideration reads 

ip[A]^Pk = Xk[A]^bk , (1) 

Q-i, where ip[A] — ilf) + g^f^ is the massless Euclidean Dirac operator. The couphng constant is denoted by g and the 
^ are the generators of the gauge group. The distribution of the color gauge fields ^ is given by the Euclidean QCD 

• • : partition function. As these gauge fields vary over the ensemble of gauge field configurations, the eigenvalues fluctuate 
about their average positions. The average spectral density is defined as 

^1 p(A) = (E^.5(A-A.[A]))^. (2) 

The average has to be performed over all gauge field configurations. 

In contrast to most other systems, however, there are two different regimes in QCD spectra which can be addressed 
in an RMT approach, the microscopic region and the bulk region. Since the Dirac operator only couples states of 
opposite chirality, the eigenvalues are pairwise positive and negative. This is the reason why two types of spectral 
fluctuations can be distinguished, namely spectral fluctuations in the microscopic limit near zero virtuality, A = 0, 
and in the bulk of the spectrum. 

Concerning the microscopic region, chiral Random Matrix Theory (chRMT) |^ incorporates the global symmetry 
properties, in particular chiral symmetry, oi ip. It predicts level repulsion between positive and negative eigenvalues 
which results in a distinct behavior of the eigenvalue density and correlations near the origin. It is possible to calculate 
spectral correlators analytically in the microscopic limit 0-0] , and to compare the predictions of chRMT with comnlete 



spectra of the lattice QCD Dirac operator on reasonably large lattices. Indeed, remarkable agreement is found |6|-|1C|] 
at the edge of the spectrum. 

Sufficiently far away from the origin, however, the repulsion of negative and positive eigenvalues should become 
unimportant. Therefore the chiral structure of the theory is not expected to be of relevance in the bulk of the 
spectrum. This is the region we will address on in this work. By comparing with lattice data, it has already been 
shown that conventional RMT properly models these fluctuations in the bulk . It is important to go beyond 

these statistical analyses made so far in order to see to what scales RMT does apply. The identiflcation of such scales 
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gives a fundamental insight into a system. Investigations of this type have been performed in great detail in disordered 
systems and in many-body systems. There, the Thouless energy Ec or the spreading width F determine the scale 
Ec/D where D is the level spacing, in which the fluctuation are of RMT type Beyond this scale, deviations 

from the RMT behavior occur, see Ref. for a detailed discussion and further references. 

Recently, theoretical studies JT5|~p7[ established a link between disordered systems and QCD. The range of validity 
of RMT, Armt, was introduced as an equivalent of a Thouless energy Ec. A scaling \rmt/D cx \/V was proposed 
where V is the four-volume of the system. Indeed, such a scaling behavior was found very recently in the microscopic 
region for deviation from RMT behavior. As argued in a corresponding effect should also be seen in the bulk 
of the spectrum. This is what we will investigate. 

The identification of this scale Armt / D in QCD spectra could lead to an improved understanding of certain features 
of QCD and allows us to separate the stochastic noise of the short range fluctuations from the true dynamics of the 
QCD vacuum. Eventually, it could be possible to set up effective models or simplify the presently used simulation 
algorithms in lattice gauge theories. In this work, we search for such a universal scale Xrmt/D in the bulk region by 
analyzing lattice data. In contrast to the microscopic region, the bulk of the spectrum is expected to have a translation 
invariant analogue of the Thouless energy. We emphasize that our analysis is self-consistent. Advantageously, it does 
not depend on any model that aims at an explanation for the occurrence of this universal scale. 

The high amount and quality of the data sets which exceed the existing ones by far enable us to considerably extend 
the energy range for our analysis. Moreover, the wealth of data makes it possible to directly address bare correlation 
functions which cannot be analyzed in most systems. Furthermore, in doing so we discuss some technical aspects, 
which are of general interest for the investigation of spectral fluctuations. 

This paper is organized as follows: In Sec. |ll| the data under investigation is presented. A detailed analysis of 



the statistical properties is given in Sec. III. This includes the introduction of the numerical unfolding approaches, 
a statistical analysis of the nearest neighbor spacing distribution, two-point spectral correlations and higher order 
spectral correlations. Deviations from the RMT predictions are found and interpreted. Summary and discussion are 



given in Sec. IV 



II. COMPLETE DIRAC SPECTRA IN SU(2) GAUGE THEORY 

The computation of large ensembles of complete spectra of the Euclidean Dirac operator for staggered fermions 
in SU(2) gauge theory has recently been performed Ref. expanding the numerical work of ||2^. In lattice gauge 
theory simulations one generates a sequence of gauge field configurations distributed according to the Boltzmann 
weight. On each of the gauge field configurations the eigenvalue equation (|^) is solved numerically on the lattice and 
a distinct partition of eigenvalues is obtained. The lattices have the size V = L"^ where L denotes the length of the 
Euclidean box with a lattice spacing a. Parameters and statistics of the simulation are summarized in Table^ The 
choice of SU(2) as the gauge group implies that every eigenvalue of ip is twofold degenerate due to a global charge 
conjugation symmetry. The random-matrix ensemble for this situation has symplectic symmetry and is referred to 
as chiral Gaussian Symplectic Ensemble (chGSE) 1^,01 . In addition, the squared Dirac operator —p^ couples only 
even to even and odd to odd lattice sites, respectively. Thus, —p^ has V/2 distinct eigenvalues. We use the CuUum- 
Willoughby version of the Lanczos algorithm to compute the complete eigenvalue spectrum of the sparse hermitian 
matrix —p^ in order to avoid numerical uncertainties for the low-lying eigenvalues. There exists an analytical sum 
rule, tr(— = AV, for the distinct eigenvalues of —p'^ We have checked that this sum rule is satisfied by our 
data, the largest relative deviation was about 10~^. 

Examples of the spectra are shown in Fig. |l|, where the average level density and the integrated average level density, 
see Eq. (|^), for a 16**, i.e. L = 16a, lattice are shown. It should be pointed out that due to the V/2 = 32768 distinct 
eigenvalues of each configuration there are millions of eigenvalues at our disposal. We used two different values of the 
gauge coupling /3 = A/g'^ where the weak coupling regime of SU(2) sets in and where most of the scaling test have 
been performed so far. Finally, the chiral condensate was obtained by fitting the spectral density and extracting p(0). 
Our findings [p^ arc in rough agreement with the values obtained by Hands and Teper [p3| for the same simulation 
parameters in SU(2) but only the 20 smallest eigenvalues have been computed by these authors. 



III. DATA ANALYSIS 



In this section we give a detailed analysis of the data introduced in the previous section in the b ulk o f the spectrum. 
We start with a description of the numerical unfolding approaches and their properties in Sec. [II A . After a short 
discussion of the nearest neighbor distribution in Sec. IIIB, we present data for spectral two-point correlations at 
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large scales in Sec. IIIC, From this we identify the equivalent of a T houless energy. Furthermore we discuss higher 
order correlators, in particular three-point correlations in Sec. Ill D . A q ualita tive explanation of deviations from 
RMT predictions which are not due to the Thouless energy is given in Sec. [HE. 



A. Unfolding 

As RMT is capable of making predictions for the fluctuations on the scale of the mean level spacing, one has to 
remove the influence of the level density by unfolding the spectra. The cumulative spectral function 

.A VII 

N{X) = / dA'^(5(A'-A0 , (3) 

J-oo 

is the number of levels below or at the energy A. It is frequently referred to as staircase function. It can be separated 
into an average part A^ave(A), whose derivative is the level density, and a fluctuating part Affluc(A), 

iV(A) = A^avc(A)+iVfiuc(A) . (4) 

The average part is determined by gross features of the system and has to be removed. The fluctuating part is in 
all relevant systems of order 0(1) and contains the correlations to be analyzed. After extraction of the average part 
-^ave(A), it is unfolded from the spectra by the introduction of a dimensionless energy variable 

^z=^avc(A.). (5) 

In this variable, the spectra have mean level spacing unity everywhere, 

l/pavc(C) = 1, (6) 

where Pavc(C) — dNg^^ciO/d^- However, the extraction of -/Vavo(A) from the data is non-trivial in our case because 
little is known analytically about the level density of QCD spectra, particularly in lattice calculations. We thus have 
to resort to phenomenological unfolding procedures. Faulty unfo lding leads to w ron g resul ts, especially on such large 



energy scales that we are interested in. In the subsequent Sees. Ill A 1 , [II A 2 and [II A3, we discuss three different 



procedures used here, ensemble unfolding, configuration unfolding and windowing, respectively. 



1. Ensemble unfolding 



In RMT one deals with an ensemble of matrices, where the matrix elements of each member are chosen randomly. 
Spectral observabl es pr edicted by RMT are calculated as an average over the ensemble. This ensemble average is 
denoted by a bar (...). But observables can also be calculated as spectral average, i.e. one performs a running 
average over overlapping intervals [Q;,a + L] of length L in the spectrum of one member. In order to distinguish it 
from ensemble average, we denote spectral averaging by angular brackets (...). In the limit of large matrix dimension 
both averages are equivalent This property is called ergodicity. 

In most experiments, one measures one - preferably long - spectrum. Thus observables are usually calculated from 
spectral average. One uses the theoretical concept of ergodicity to compare the RMT predictions with the experimental 
results. In our case, however, the data consists of configurations, i.e. forms an ensemble. Hence, questions related 
to ergodicity arise not only for the calculation of observables, but also in the determination of the staircase function, 
i.e. in the unfolding procedure. We have in principle two very different ways of unfolding our data: first, ensemble 
unfolding, -/Vavc(A) = N{X), i.e. we determine the smooth part of the staircase function by averaging over the ensemble, 
and second, configuration unfolding, A'avc(A) — {N{X)), i.e. we determine the smooth part of the staircase function for 
every configuration separately. The results differ considerably, iV(A) ^ {N{X)), for most of the configurations. The 
ensemble averaged staircase A^(A) for the lattice QCD Dirac operator is shown in Fig. |l|. We find iV(A) by dividing 
the energy range in m bins with width AA and average the density p(A, A + AA) for each bin over all configurations. 
We then calculate the staircase function as N{X) = Pi^i' + ^A)AA, with Am = A. In Fig. ||, the difference 

between the ensemble averaged staircase function and the configuration wise averaged ones, for V — 16"* and /3 = 2.4, 
for 50 arbitrarily chosen configurations is plotted. Each data point represents the difference N{Xij) — {N{Xij)), where 
i enumerates the eigenvalues, i = 1, . . . , 32768 and j is the configuration number j = 1, . . . , 50. We plot only every 
500th eigenvalue. There are deviations of about N{X) — {N{X)) = 0(10^) in certain energy ranges. 
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If the spectra are unfolded using the ensemble averaged staircase function N{X), observables should then also be 
calculated as an ensemble average for a fixed value of A. But we checked that our results do not depend on A in a wide 
range of the bulk. This property is called translational invariance. It is actually not present in the microscopic region, 
where it is destroyed by point wise symmetry between positive and negative eigenvalues fl^ . Translational invariance 
in the bulk allows us to calculate observables from running average over overlapping intervals for each configuration. 
We choose an overlap of 90% for two consecutive intervals. Then we average over all configurations. This improves 
the statistics of the result considerably. 



2. Configuration unfolding 

We now unfold each configuration separately. Observables are then calculated for each configuration by running 
spectral average. Thereafter we average over the ensemble. The basic characteristics are already obtained for one 
single configuration, though the statistics is considerably improved by ensemble averaging. This is in the same spirit 
as it was done in spectra of nuclei and complex atoms [p5| . These spectra were unfolded for each nuclei or atom 
separately. Then observables were calculated as described above, i.e. first taking the spectral and then the ensemble 
average. In this case the ensembles consist of nuclei or atoms of different types. 

Configuration unfolding is, in contrast to ensemble unfolding, not a unique procedure. One has to find either fits 
to the average staircase function or to smooth it in some way. A priori, there is no criterion whether the numerical 
estimated iVavc(A) is close to the real one or not. Thus, we use three different approaches and carefully compare them 
with one another to eliminate as many sources for mistakes as possible and to obtain consistent results. 

First, we fit A^(A) to a polynomial of degree n, 

n 

3=0 

where n is a small integer, n = 2, ... ,5. This approach is motivated by the fact that almost all physical systems 
are known to have a level density which is as smooth as a polynomial. In our case this ansatz is supported by 
pertubative calculations. Strong coupling expansions for SU(2) with staggered fermions have been performed |2q ] 
and, furthermore, 1/Nc expansion of the QCD level density p7| , both suggesting a smooth level density. The former 
gives a semi-circle whereas the latter explicitly predicts a polynomial increase. 

Second, we use the Gaussian method which was originally developed by Strutinsky ||2^ . One replaces the (5-functions 
in Eq. {m by Gaussian functions with a width A which yields a smoothed staircase 



/\ max 

The summation runs from the smallest eigenvalue X^nin to the largest eigenvalue Amax in the interval under consid- 
eration. The limit A ^ restores ^-distributed eigenvalues, whereas the fiuctuations are smeared out for finite A. 
The optimal parameter Aopt is found by a x^~fit of Na{X) to N{X). Then we identify 

(iV(A)) =iVA„p.(A) . (9) 

Third, we perform a local unfolding by calculating the unfolded eigenvalues directly with the formula 

- 6 = , (10) 

with local mean level spacing 



i+k 



j=i-k 

Here 2k is the number of consecutive level spacings over the running average is performed. 

Whatever approach one decides to use, a necessary condition is that on the unfolded scale the average number of 
levels in an interval of length L should equal this length. This is a very important requirement because we are also 
interested in very large energy scales. This assures that the spectrum on the numerically constructed dimensionless 
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scale 5 has mean level density unity. Consider the interval + L] which contains na{L) eigenvalues. Spectral 
average (. . .) and ensemble average (. . .) have to yield 

ML))=L. (12) 



In Fig. |3| the difference between the calculated mean number of eigenvalues {na{L)) and L is plotted as a function 
of L. For the Gaussian method the difference appears to be zero to all scales. While it is small and does only appear 
at large L for polynomials fits with n> i, stron g deviat ions from the zero line already appear at small L for n — 2. 
In the case of local unfolding, the difference L — {na{L)) is positive for small k, i.e. there are on average less levels in 
a given interval than there should be. For growing fc, it becomes negative with ever stronger deviations from the flat 
line. The averaging interval length k for which the difference equals zero is fc « 100 for V — 16'*. For other lattice 
sizes this averaging interval is slightly smaller. We take this as the optimal parameter for this approach. From Fig. ^ 
we learn, that the necessary condition (n2|) is fulfilled only for the Gaussian approach, polynomial fit with n > 3 and 
local unfolding with k « 100 for V = 16 , all other choices of the parameters must be rejected. 

It should be mentioned that a new artificial scale both for local and Gaussian approach is introduced, namely the 
averaging interval length k and the width A, respectively. Therefore, one should be cautious in the interpretation of 
effects seen on scales L larger than current value of the corresponding parameter. In units of the mean level spacing 
D we find a width of the Gaussian as A/D « 100 at a 16^ lattice. On the other hand, both approaches have the 
advantage that no particular function for the average level density has to be assumed. 

We checked all our numerical unfolding approaches with the spectrum of a very different system. We used the 
spectrum of quantum chaotic billiard that was simulated in a microwave experiment. In billiards, the Weyl formula 
gives an analytically expression for the mean level density p9[ . With our phenomenological approaches, we indeed 
obtained quantitatively the same results. 



3. Windowing 

Ideally, an unfolding procedure should only remove the global variations of the spectral density, i.e. in our case the 
overall behavior seen in Fig. |l|. For reasons which will become clearer later, it is difficult to numerically distinguish 
the global variations from the local fluctuations. This is in particular the case for data of large lattice size. In other 
words, we might have removed too much by some of the unfolding procedures, while we might have removed too little 



by others. This will be discussed in great detail below, especially in Sec. IIIE. 

One has to ensure that any deviations seen in the spectral statistics arc not due to global variations in the average 
density which were not removed adequately. One way is to use different unfolding approaches and compare the results 
carefully. Another way is to take only a small window of the spectra in which the global variation of the density is 
expected to be small. Thus we choose an interval [A, A + SX] and calculate the ensemble averaged mean level spacing 
D for it. Wc then rescale the eigenvalues in this interval as 

^,^XJD , A<A, <A + ^A. (13) 

This is done in the same manner as in Q where the microscopic region was considered. However, it is not clear 
beforehand that a scale in the spectra, if any, does not exceed the interval length SX. Unfolding, if done correctly, 
allows to make investigations to much larger scales. 

This approach is closely related to ensemble unfolding defined above. Indeed, the results coincide, but with a less 
statistical significance by only rescaling. By using this approach, we intended to avoid any unfolding procedures. As 
we will see later, there are a slight, but still systematic variations of the spectral density within the small window. 



B. Nearest Neighbor Spacing Distribution 

The nearest neighbor spacing distribution P(s) probes the fluctuations on short scales in the spectra. It is the 
probability of flnding the distance s between adjacent level on the unfolded scale. It contains all correlations of order 
k > 2. In the case of completely uncorrelated levels which is referred to as Poisson regularity ||l|], it is given by 
P{s) = cxp(— s). In the case of GSE type correlations, Wigner surmized the shape 
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which is very close to the exact GSE result. As shown in Fig. ^ the data is in perfect agreement with the prediction. 
This is as well true for the large lattice, V = 16^ (left part), as for the small lattice, V = A"^ (right part). The 
spacing distribution of the intermediate lattice sizes are not distinguishable from the both shown. We have complete 
agreement between theory and lattice data for any choice of lattice size and coupling constant. 



C. Two— Point Spectral Correlations 

In an interval of length L in units of the mean level spacing, the mean number of eigenvalues should be equal to L, 
see Eq. (|l|). The variance of this number is defined by 



I]^(L) = ((L-n4L))2) . (15) 

Thus, an interval of length L contains on average L ± ^JT?{L) levels. For uncorrelated Poisson spectra S^(L) = L. 
RMT predicts for the number variance stronger correlations, namely Y?{L) ~ logL. 

Another important quantity is the spectral rigidity A3(L), defined as the least square deviation of the staircase 
function from the straight line ||3^ , 

A3(L) = ^^min^^B ^"^"^ d^{N{^) -A£,- Bf^ . (16) 

Since it can be expressed as an integral over the number variance 

A3(i) = dr{L' - ^L'r + r')j:\r) , (17) 

it is smoother than E^(L). 

The number variance can be expressed as an integral of the two-point cluster function Y2{r) [ ^0| , which depends 
for translational invariant spectra only on the difference "r = |^2 ^ Ci| between two levels at and ^2, 

J:^{L)^L-2[ {L-r)Y2{r)dr . (18) 
Jo 

The cluster function is related to the two-point correlation function X2 (r) which measures the probability density to 
find two levels at a distant r by X2{r) = 1 — 12(7'). In contrast to P(s), these two levels are not restricted to adjacent 
ones. 

In Fig. ^ number variance and spectral rigidity for scales up to L = 20 and L = 100, respectively, are shown. 
Lattice data and RMT predictions agree remarkably well, even the oscillations in Y,2{L) are accurately reproduced. 
Naturally, previous analyses ||ll| with smaller data sets have less statistical significance. Two-point cluster and 
correlation function which usually are not accessible in data analysis are shown in Fig. ^. Again, the agreement is 
impressive. 

Beyond this scale there are considerable deviations of /S.j,{L) as well as of S^(L) from RMT predictions which 
depend on the unfolding procedure used. We mention that on general grounds one can show that any scales in Y?{L) 
and A3(L), say Lp and L^, respectively, are related by : = 4 : 1, or so see Fig. ^ Thus, any deviations 
from RMT behavior appear at smaller L in S^(i) as compared to A3(L). 

If we unfold with the ensemble staircase function N{X), we obtain the following results. The number variance can 
be seen in Fig. |^. Data for different lattice sizes and different gauge couplings are shown in comparison to the RMT 
predictions for some regions of the spectra. We find that the point where the deviation sets in, scales with the square 
root of the lattice volume, 

(19) 

The numerical constant C is approximately given by C « 0.3. This should be compared with the result obtained in 
p8[ for the microscopic region. There, the scaling Armt/-D ^ 0.3 .. . Q.7\/V was found. This is independent of the 
region of the spectra we consider and of the coupling strength (3. This is shown in Fig. 0. There different regions 
of the spectra are considered, each corresponding to different values of p(A), see Fig. |l|. The results are the same. 
Furthermore, the deviation points for different (3 appear to coincide, whereas the local average density depends on 
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the gauge coupling, p{\) = p(A,/3), see Fig. |[ The scahng relation ( p^ can nicely be seen from Fig. ||, in which the 
L axes of Fig. ^ is rescaled with VV. We see that the crossover from RMT to non-universal behavior appears to be 
the same for all lattice sizes independent of /3. But the slope varies for different couplings and regions of the spectra. 
When we use windowing instead, we get the same results as obtained by ensemble unfolding. But the data points 
scatter more compared to Figs. ^ and g. 

The importance of a proper choice of the unfolding method becomes manifest as the above picture changes drastically 
if we unfold each configuration separately. As displayed in Fig. |^, the polynomial unfolding leads to an overshooting 
of the data over the predictions but further out, compared to ensemble unfolding, while in the Gaussian as well as in 
the local case saturates. Note the different scale compared to Fig. |^ and also the difference in scale between 

number variance and spectral rigidity, as mentioned above. The result of the polynomial approach does not depend on 
the degree n of the polynomial. Furthermore we find no scaling with \/V for the deviations of polynomial unfolding, 
see Fig. |l0|. The deviation point appears to be the same for different lattice sizes. The saturation of the small lattices 
is due to the limited number of eigenvalues in the considered energy range. The same picture arises for the number 
variance: overshooting for the polynomial, saturation for Gaussian and local approach. The general tendency of this 
results are already obtained for each configuration separately, but the data points scatter. After averaging over all 
configurations the scattering becomes much smaller, see Figs. |^ and ^. 



D. Higher Order Spectral Correlations 

The wealth of data allows us to go beyond a previous analysis of higher moments of the eigenvalue partition, 

(20) 



We notice that /^2(-^) — The skewness and the excess PQ] are defined by 



and 



7i(L)=/X3(L)/X2(L)-=^/' 



3 , 



(21) 



(22) 



respectively. The comparison of RMT predictions with lattice data for these both quantities in Fig. 11 again shows 
very good agreement. 

The measures ji{L) and 72 (i) only contain a small amount of information of the spectral correlations. Moreover, 
71 (L) and 72 (i) also involve lower order correlations: 71 is a combination of the two- and the three-point correlator, 
Y2(r) and Y3{ri^r2), and 72 involves in addition the four-point correlator i4(ri, r2, rs). The representation of the 
moments in terms of integrals over the correlators reads 



H3{L)=L-6 {L - r)Y2{r)dr - 6 dri 



dr2{L - ri - r2)Y3{ri,r2) 



(23) 



and 



fi4{L) ^ L - {U ~ 12L) [ {L - r)Y2{r)dr + 12 [ 
Jo Jo 



{L-r)Y2{r)dr 



+ 36 
- 24 



dri 



dri 



dr2{L - ri - r2)F3(?'i, ''2) 



dr2 



drsiL - ri - r2 - r3)y4(ri, r2, rs) 



(24) 



Obviously, by analyzing 71 (i) and 72 (L), one cannot easily estimate to what extent the three- and the four-point 
correlators themselves obtained from the lattice calculations follow the predictions of RMT. 

Here, the three- and the four-point cluster functions, Y3{ri,r2) and 14(^1, 7^2, ra) are written as functions of the 
arguments (i = 1, 2, 3) which are defined terms of the original arguments (i — 1, 2, 3, 4) by 



ri^^2-^l, ^-2=6-6, ^'3=^4- 6 



(25) 
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We constructed from the data the two-point and the three-point correlation functions X2{s) and X^{si^S2) and 
the corresponding cluster functions Y2{s) and 13(51,52)- Here, for convenience, we redefined the arguments for the 
three-point correlators as follows: 

52 = ^3 - Cl = ^1 + ?'2- (26) 

The results for ^3(51, 52) and 13(51, 52) are plotted with error bars in Fig. |l2| for some given values of 5i. The results 
do not depend on unfolding. In the construction of these correlators, we first performed a spectral average by using 
the translational invariance due to unfolding, and then averaged over the ensemble. The errors for ^3(51,52) and 
1^5(51,52) were estimated as the variance of the ensemble fluctuations. Once more, very good agreement with the 
RMT predictions for ^3(51, 52) is found, apart from a small systematic deviation which we believe can be understood 
as follows. From the relation between ^3(51, 52) and ^2(5) 

^3(51, 52) = ^3(51, 52) - 2 + ^2(51) + ^2(53) + X2(5i - 52), (27) 

one has 

X3(5l,52)|s,-.oo =^2(51). (28) 

Therefore, even a small point -deviation of the two-point correlator at 51 from the theoretical predictions can result 
in a systematic deviation of the whole curve ^3(51, 52) versus 52 for this given 5i. For ¥3(51, 52), the quality of the 
agreement with RMT is only slightly reduced, but still remarkable. In addition to the systematic deviation, one can see 
some random fluctuations around the theoretical curves. This is because ^3(51, 52) is the disconnected correlator, and 
it should represent the true three-point correlation, while ^3(51, 52) also contains the two- and one-point functions. 
They play a dominate role and are in good agreement with the corresponding GSE predictions. We notice that this 
analysis was only possible due to the extremely higher number of levels available. 



E. Qualitative discussion of the deviations for configuration unfolding 



Using ensemble unfolding, we find deviations from RMT behavior, which scale with the square root of the volume 
according to the theoretical predictions. But as we unfold each configuration separately this effect vanishes. There 
are still deviations left but none of them show a ^/V scaling law. Moreover, we have a dependence of the results on 
the unfolding approach. 

Concerning local and Gaussian unfolding, an explanation seems to be easy at hand. As mentioned above, both 
procedures have an intrinsic, artificial, scale. In units of the mean level spacing it has the value L « 100 for V = 16^ 
in both cases. This is approximately where the saturation of the statistics seen in Fig. ^ sets in. We conclude that 
this artificial scale causes S^(L) and A3(L) to saturate. In other words, both approaches are not capable to allow 
statements at scales L ^ 100. Nevertheless, both approaches do not show a behavior as shown in Figs. |^,^ for L < 100. 

On the other hand, the polynomial unfolding also deviates from RMT predictions, see Fig. |l^. As this approach 
does not contain an additional scale, we can rule out effects like the one discussed above. The question is, whether 



the deviations in Fig. 10 are due to a Thouless energy or if they have another origin. 

The fluctuating part of the integrated level density Affluc(A) should be of order 0(1), as mentioned above. In the 
upper part of Fig. O the difference between the real staircase function and the smooth polynomial staircase function, 
N(X) — A'poiy(A), for one specific configuration is plotted. This picture remains qualitatively the same whatever 
configuration is chosen. The polynomial fits have a systematic deviation from a smooth behavior larger than 0(1). 
The difference between the ensemble staircase and polynomial fit N{\) — A'^poiy (A) is shown in the lower part of Fig. HS. 
A polynomial of degree n — 3 gives the same result as n = 4. To obtain a better insight in this obviously not universal 
behavior, we calculate the power spectrum 

F{f) = J dAe2-/^X(Ai, A2, A)— (iV(A) ~ iVpoiy(A)) . (29) 

By construction, the derivative in the integrand gives the fluctuations of the level density. The window function 
^(Ai, A2, A) has to be introduced because we only have a finite interval of eigenvalues in the Fourier transform over 
the whole real axis. It is zero outside the interval Ai < A < A2. The choice is not unique inside |^l|. Thus, the Fourier 
transform is a convolution of the transforms of N{X) — 7Vpoiy(A) and i^(Ai, A2, A). In order to reduce the influence of 
the Fourier transform of i4r(Ai, A2, A) on the results as far as possible, we use a triangle window [Q. The result is 
shown in Fig. |lj. In the right part one sees the Fourier transform of the ensemble averaged density. Only the very first 
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peak is left, both for polynomial of degree n = 4 and n = 5. The latter is reduced in amplitude. On the left side the 
transform of an arbitrary chosen configuration is plotted. The first peaks are reduced in amplitude again for n = 5, 
whereas the remaining ones are the same for both degrees. We conclude that only the first peak, corresponding to the 
long wave part of Fig. is common to all configurations. All others fluctuate from configuration to configuration. 

We conclude that the average level density Pavc(A) and thus the average integrated level density N^vcW consist of 
two parts, namely a very smooth polynomial-like part and another, non-universal, part, 

A^avc (A) = Npoly (A) + TVosc (A) . (30) 

We stress again that the existence of a polynomial-like smooth part is suggested by pertubative and 1 /Nc expansions 
of the QCD level density pq,M. It is expected to be governed by the available phase space: for free fermions the 
spectral density is given atfAp— > oo by Pavo(A) = NcV\X\'^ / in'^ [Q, where Nc is the number of colors. This also 
holds in a 1 /Nc expansion of the interacting theory on scales which are large compared to the hadronic scale . 
This is the region we investigated in the spectra, as the eigenvalues are given in units of the inverse lattice spacing 
as ~ 10 fm""'^ « 2 GeV, see Figs. |l| and This is why we tried to approximate the average level density by a 
function which is as smooth as a polynomial. 

The additional structure of the level density appears to have similarities to oscillations, see Figs. |l3| and Thus, 
we refer to it as "oscillatory part", 7Vosc(A). This oscillatory part explains the different behavior of A3(L) and S^(L) 
for large L for different unfolding methods. The polynomial unfolding is clearly unable to remove the oscillations 
and fits only iVpoiy(A). Thus, the oscillatory part is still present in the unfolded spectrum. The presence of these 
oscillations leads to values for A3(L) larger than predicted by RMT, because a fit to a straight line can only be done 
in a less accurate manner. In contrast to that, the Gaussian and the local unfolding is capable to fit iVpoiy(A) and 
part of A'osc(A). However, it is not clear whether the fit to the oscillatory part is done completely or if, on the other 
hand, it does not smooth out part of the universal fluctuations, i.e. overfits the data points. But, because of the 
saturation of S^(i) and A3(L), see Fig. ^ we think that probably the latter happens. 

However, as argued above, since we expect the physical density to be as smooth as a polynomial, the oscillatory 



part is likely to be a lattice artifact. This is suggested by Figs. 13 and ^ which show that these oscillations live 
on the scale of the inverse lattice spacing 1/a. As |A| cannot be arbitrarily large in lattice gauge theory, due to an 
ultraviolet cutoff in momentum for finite lattice spacing a, the increase of the density is disturbed by lattice artifacts. 
For the large lattices, i.e. V = 16^, the deviations due to lattice artifacts set in at approximately the same scale as 
the expected equivalent of a Thouless energy. This can be seen from the ^/V scaling of the smaller lattices. A rough 
estimate gives Xrmt/D « 30 for V — 16^ for S]^(L). Therefore one should be careful with the determination of the 
deviation point for large lattices in the bulk of the spectra. After removing this part from the data, we find for (L) 
reasonable agreement with the RMT prediction up to at least L « 300 and with some uncertainties even to L « 500, 
see Fig. ^ Because of the complicated structure of the average level density it is not possible to make any statement 
beyond this scale. In chaotic billiards, so-called bouncing ball modes generate effects which are similar to the ones 
here [p9| . In that case, however, an analytical prediction for the oscillatory behavior was at hand. Such a result is 
also highly desirable in our case. It would be needed to furnish our phenomenological removal of the oscillations with 
a theoretical justification. 

In any case, we may use the information displayed in Fig. |l^ to phenomenologically remove the oscillatory part. We 
cut the power spectrum at a certain interval, back transform the remaining peaks in the frequency interval [0, /cut] 
and subtract the smooth oscillatory part of the integrated level density obtained in this manner. We do this for 
each configuration separately. This procedure changes the result of the statistical analysis in a crucial way, as can 
be seen from Fig. There, the results for spectral measures are shown for different choices of /cut- Because only 
the very first peak is common to all spectra whereas higher frequencies appear to be specific for each configuration, 
a choice of /cut > 3.0 • (2a) is actually too large. This manifests in a saturation of the statistics in Fig. ^ for 
/cut = 7.0 • (2a), 10.0 • (2a). A cutoff of /cut = 1.5 . . .3.0 • (2a) seems to be the best choice, but we are not able to 
give an exact value. This figure also shows a comparison between /3 — 2.4 and /3 = 2.5 for V = 16^. Both values of 
(3 give almost the same result. The procedure also removes the slight differences seen in Fig. ^ for the polynomials 



for L > 200. From all this, we conclude that the deviations in Fig. 10 are due to a non polynomial-like part in the 
average level density and not due to an equivalent a Thouless energy. 



A possible way to circumvent the problems encountered by lattice artifacts is windowing as discussed in Sec. HI A 3 . 
However, this works only if the size of the windows SX is at least so small that the oscillations of Fig. |l^ are not 
important, i.e. SX ^ a^^. Another way might be ensemble unfolding, Sec. Ill A 1 . As the lattice artifacts are seen in 



each configuration as well as in the average integrated level density this approach might remove the oscillations. But 
it is by no means clear if it really does, especially for large interval lengths. As both approaches give similar results for 
small lattices, V < 10*, for intervals L ^ 20, we conclude that the artifacts are not important on these scales. This is 
supported by the observations that deviations due to them set in at L w 25 for polynomial unfolding, independently 
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of the lattice size, see Fig. |l^. Thus, these artifacts become important for analysis of larger lattices, i.e. V = 16^, and 
for even larger ones in future examinations. 

From all this we conclude that we do not see an equivalent of a Thouless energy if we unfold each configuration 
separately. This is in complete contrast to the results gained by ensemble unfolding. The former surely removes the 
fluctuations in Fig. ||, whereas the latter does not. We conclude that this fluctuation with respect to the ensemble 
already contains the information needed to determine the Thouless energy. After removing this information it seems 
that there is no further information in the spectra, i.e. we flnd agreement with RMT on huge scales. 

IV. SUMMARY AND DISCUSSION 

After summarizing our results in Sec. IV A, we discuss our flndings in Sec. [V B| . 



A. Summary 

We presented a detailed analysis of statistical properties of complete eigenvalue spectra for staggered fcrmions for 
SU(2) lattice gauge theory for various couplings and lattice volumes. Unfolding the data posed certain difficulties which 
are not encountered in other systems. The staircase function found by an average over the ensemble differs in most 
cases from the smooth staircase of one specific configuration. The deviations are as large as (7V(A)) — N{X) = 0(10^) 
for V — 16''. Varying the unfolding approach leads to different results for large scales. In particular, there is a drastic 
difference between ensemble and configuration unfolding. 

Using ensemble unfolding, we find a range of validity of RMT, giving \rmt/D = CW . The numerical constant 
is approximately given as C ~ 0.3. which is compatible with the result obtained in for the microscopic region 
where the scaling X-rmt/ D ~ 0.3 . . . Q.7^/V was found. The same results are obtained if we use windowing, but with 
less statistical significance. 

By unfolding each configuration separately, we do not see any scaling of this type. This procedure obviously removes 
the fluctuations of the staircase function relative to the ensemble average seen in Fig. ^ We notice that fluctuations 
over the ensemble can already be observed in the integrated level density. But there are still deviations depending on 
the unfolding approach used. In particular there is an overshoot for polynomial unfolding at large L. The reasons for 
this lies in special features of the average level density, which consists of at least two parts: a smooth polynomial-like 
and an oscillatory part. Hence the Thouless energy is due to the fluctuations in the ensemble. 

We want to clarify the difference in the notion of ergodicity used in spectral analysis and in lattice calculations, 
respectively. Concerning the analysis of spectra, a system is called ergodic if spectral and ensemble average yield 
the same results. It can be proven rigorously that random matrices have this property in the limit of large matrix 
dimension. However, here we face a different situation. We emphasize that these issues do not affect the notion of 
ergodicity as used in connection with lattice simulations. The latter is defined by the requirement that the space of 
possible lattice configuration is explored reasonably fast by the simulation algorithm. We emphasize that the time 
history of the spectra is not of importance for the calculation of the spectral observable of Sec. |l|, because it is only 
relevant that the configurations used, see Table |, are representative members of the ensemble, i.e. ergodicity in the 
sense of lattice simulations is fulfilled. This is ensured by our algorithms as described in Sec. ||. Furthermore, this is 
supported by the observation that our results do not change if we take arbitrary subsets of configurations, only the 
scattering of the data points increases. 

Furthermore we analyzed the nearest neighbor-spacing distribution, skewness and excess for small L. Due to the 
large amount of data, our results achieved a quality never reached before. Due to translational invariance in the bulk of 
the spectra the general tendency for all of the observables can already be seen for one specific configuration. Averaging 
over all configurations increases the statistical significance of the results. The results for small values L considerably 
improve the statistical significance of previous analyses |ll|,|l2|. To the best of our knowledge, we presented the first 
statistically highly significant analysis of bare two- and, importantly, three-point spectral correlators. Very good 
agreement between lattice data and RMT predictions is found. We find minor deviations for the three-point cluster 
functions. In our opinion, their origin are probably very small point-deviations of the two-point functions. 

B. Discussion 

Some scenarios for the level statistics are shown in Fig. |l^. There the spectral rigidity is plotted. The two solid 
lines represent Poisson and RMT behavior, A3(i) ~ L and A3(i) ^ logL, respectively. The dashed lines represent 
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schematically two different cases. Case 'A' correspond to a scenario known from systems with few degrees of freedom. 



The shortest periodic orbit in phase space sets a scale which forces A3(L) to saturate at a certain La l33|. But by 
increasing the degrees of freedom of the system, this scale become ever larger and is hard to be seen. As we may view 
lattice gauge theory as a many body problem or disordered system, we do not expect to see this scale in the level 
statistics. Therefore one should not be tempted to conclude that the saturation seen in Figs. is the analogue of 
"saturation of fluctuation measures" which was observed in the energy level statistics for classically chaotic systems 
and was interpreted by Berry in terms of shortest periodic orbits . Unlike in the case of quantum billiards where 
an analytic expression for the average spectral density exists, the A'avc(A) that we obtained in our analysis was only 
a numerical result. Most likely all three approaches, i.e. Gauss, local and modified polynomial, fit also parts of the 
universal fluctuations which shows up as a saturation of the statistics A3(L) and S^(L). 

Another scenario is given by case 'B' which lies between Poisson and RMT behavior. It can arise in three different 
physical situations [Q: 

1. Systems in few degrees of freedom between regularity and chaos. The spectral rigidity lies between the Poisson 
limit, which applies to regular systems and the RMT limit which applies to chaotic systems, provided the scale 
La is much larger than Lb- 

2. Disordered systems in d dimensions. The time tc for the classical diffusion through the system of size where 
L5 is the size in each dimension determines an energy scale 

Ec - n/tc , (31) 



the Thouless energy. This in turn sets the scale \bmt/D where D is the single particle mean level spacing. 
For energy separations in units of D smaller than this scale, RMT fluctuations are seen. For larger energy 
separation, deviations as sketched in Fig. ^ occur. 

3. Many body systems such as nuclei, molecules or complex atoms. In the language of condensed matter physics, 
these systems are zero-dimensional. Consider the Hamiltonian of such a system H = Hq + Hi where Ha has 
a certain property and Hi breaks this property. The influence of Hi on the statistics of the spectrum of Hq is 
measured in terms of the spreading width 

r = 27T{Hf)/Do (32) 



where Dq is the mean level spacing of Hq and (H^) the mean square matrix element of Hi. In particular, if Hq 
has Poisson statistics and Hi RMT statistics the spectral rigidity acquires the form of Fig. |l^ and F sets the 
scale Xbmt/D- 

As discussed above, situation 1. is not likely to apply. We emphasize that the Thouless energy and the spreading 
width F deflned in 2. and 3. are closely related concepts. Actually, in his original paper, Thouless ||l^ deflned Ec as 
a special kind of spreading width. The precise relation between Armt/-D and Ec has been studied in great detail in 
the framework of RMT pM. We notice that the definition in 2. which is most commonly used in condensed matter 
physics, relates Ec to dimensionality, whereas F is defined in zero dimensions. Thus, the sheer existence of a scale 
Armt/-D can imply, but does not necessarily imply that the deviation from RMT is related to a diffusion in a d 
dimensional disordered system. 

Recent theoretical studies aim at establishing p^-p7| a link between disordered systems and QCD. In these works, 
a range of validity Armt of applicability of RMT was considered. Using the Gell-Mann-Oakes-Renner relation, it 
should be determined by the pion decay constant and the chiral condensate E, 

Armt ~ — 7= ■ (33) 



In the microscopic region the mean level spacing is D ~ tt/CEV), so that Eq.(g3|) can be rewritten as 

Armt 



D ^34) 

The predicted scaling behavior was found very recently in the microscopic regio n [|l8i . There a crossover from RMT 
to non-universal behavior was found with the scaling Eq. ( p^ . As argued in Rcf. |l6|, a similar relation should also be 
seen in the bulk of the spectrum. Here, however, the mean level spacing D is that of the bulk. In the present work, we 
have shown that, first, the scale Armt/-D exists in the bulk, and, second, that it shows the predicted scaling behavior. 
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The latter observation is a support for the ideas that Hnk QCD to disordered systems and thus to a diffusion in a d 
dimensions. However, as outUned above, our findings are necessary but not sufBcient for this conclusion. They do not 
rule out other physical mechanisms that lead to a spreading width or Thouless energy. This underlines the strength of 
our analysis: it does not depend in any way on a model for this mechanism. It is a self-consistent statistical method 
to find the scale Armt/-D. 

We hope that the identification of this scale Armt /D can help to improve our understanding of certain features of 
QCD as it may allow one to separate the stochastic noise of the short range fluctuations from the true dynamics of 
the QCD vacuum. Most desirably, these results could help to design effective models or to simplify the presently used 
simulation algorithms in lattice gauge theories. 
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TABLE I. Lattice parameters and statistical analysis of the complete spectra of the Dirac operator. 
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FIG. 1. Average level density p(A) for j3 = 2.4 and j3 = 2.5 (left plot) and integrated average level density N{X), see Eq.(|^, 
for {3 = 2.5 (right plot). The eigenvalues are given in units of the inverse lattice spacing (2a)~^. The bin sizes in the left plot 
is 0.01 • (2a)-\ 
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FIG. 2. Difference between the integrated level density N{X) averaged over all 921 configuration (/3 = 2.4) and real 

data. Each dot represents the value of N{\i.j) — N{\ij). Index i enumerates the eigenvalues, i = 1, . . . ,32768, and j is the 
configuration number, j = 1, . . . , 50. The 50 plotted configuration where chosen arbitrarily. Only every 500th eigenvalue is 
shown. 
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FIG. 3. Value of the quantity L — {na{L)) for the three unfolding approaches on a 16* lattice. Prom left to right the 

polynomial, Gaussian and local unfolding is shown. In the left plot diamonds are data for n = 2 and the cross and circles arc 
n = 3 and n = 4. In the right plot the data points from top to bottom correspond to an averaging interval of A; = 20,100,300 
and 900, respectively. 
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FIG. 5. Integrated two-point functions number variance E^(_L), spectral rigidity A3(L) for small L on a 16''-lattice. The 
solid line represent the RMT predictions and the dots the data. On this scale the presented data points do not depend on 
unfolding. Note the difference in the scale of the L axes between A3{L) and S^(L). 
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FIG. 6. The two-point correlation function X2(r)(left) and the cluster function Y2{r) (right) as a function of r, compared 
with the GSE predictions. The result is independent of the unfolding approach. 



18 




19 




20 



0.3 
0.2 
0.1 



^ 0-3 
< 0.2 
0.1 



0.3 



0.2 



0.1 



polynom 




Gauss 




local 



500 1000 1500 2000 
L 




200 400 600 800 
L 

FIG. 9. Deviations of the spectral rigidity A3(L) and number variance S^(L) from the RMT- predictions on a 16''-lattice 
for large L. From top to bottom the results for polynomial with degree n = 3, Gaussian and local unfolding with averaging 
interval length k — 100 are shown. Note the different scale on the abscissa compared to Fig. ^ 
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FIG. 10. Comparison between RMT and lattice data by unfolding each configuration separately with a polynomial. 
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FIG. 11. Integrated three-point and four-point function, skewness yi{L) and excess 72(i), respectively, as defined in 
Eqs.(pl|)-(|24|). The lattice size is V = 16'', as in Fig. ^. The solid line represent the RMT predictions and the dots the 
data. On this scale the presented data points do not depend on unfolding. 
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FIG. 12. The three-point correlation function X3(si, S2)(left) and the cluster function y3(si, S2)(right) as a function of S2 
for different values of si — 0.70,1.40 and 2.55 (from top to bottom), compared with the GSE predictions. As in Fig. ^ the 
results are independent of the unfolding approach. 



FIG. 13. Difference between the fitted polynom like staircase function and the real staircase function for one arbitrarily 
chosen configuration (upper part). The lower part shows the difference between the staircase found by ensemble averaging and 
a polynom fit to it. The degrees of the polynomials are n = 4, 5. Polynomial of degree n — 3 gives the same result as n = 4. 
The plotted interval contains approximately 16000 eigenvalues. 
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FIG. 14. Square of the Fourier transform of the oscillation shown in Fig. ^ but for d{N{\) - iVpoiy(A))/dA instead of 
N{X) - Afpoiy(A), as given by Eq.@. 




FIG. 15. Number variance S^(L) and spectral rigidity A3(L) for polynomial unfolding with the withdraw of the long wave 
length oscillations, as explained in the text. In each plot the data corresponds from top to bottom to a cut of /cut =0, 1.5 ■ (2a) , 
3.0 ■ (2a), 7.0 ■ (2a) and 10.0 • (2a), respectively. 
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FIG. 16. Possible scenarios for the spectral rigidity. Plotted are the linear Poisson behavior and the logarithmic increase of 
the Gaussian ensembles (GE) as predicted by RMT, respectively. 'A' correspond to a saturation due to shortest periodic orbits 
and 'B' to linear increase due the scale set by a Thouless energy. 
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This figure "diff3.gif" is available in "gif" format from: 
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